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Abstract We present a reduced basis method for the simulation of American option 
pricing. To tackle this model numerically, we formulate the problem in terms of a 
time dependent variational inequality. Characteristic ingredients are a POD-greedy 
and an angle-greedy procedure for the construction of the primal and dual reduced 
spaces. Numerical examples are provided, illustrating the approximation quality and 
convergence of our approach. 



1 Introduction 



We co nsider the problem of American option pricing and refer to (lAchdou and Pironneau , 
2005) and the references therein for an introduction into computational methods for 
option pricing. While European options can be modelled by a parabolic partial dif- 
ferential equation, America n options result in additional inequality constraints. We 
refer to dHager et all 1201 Oh for a possible numerical treatment by primal-dual fi- 
nite elements and to (iGlowinskii 12008c iGeiger and Kanzowl 120021) for an abstract 
framework on the theory of constrained variational problems. We are interested in 
providing a fast numerical algorithm to solve accurately the variational inequality 
system of an American put option for a large variety of different parameter values 
such as interest rate, dividend, strike prize and volatility. Reduced basis (RB) meth- 
ods are an approp riate means for standard parametrized parabo l ic partial differen 
tial equations, cf. (lHaasdonk and Ohlber ger. 2008; Rozza, 2005: IVeroy et a i l2003l 
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Buffa et all 1201 II) and the references therein. These are based on low-dimensional 



approximation spaces, that are constructed by greedy procedures . Convergence be- 
havio r of these procedures are known in some cases (Buffa et al, 201 1; Haasdonk, 



2011). The computational advantage of RB-methods over standard discretization 



methods is obtained by its possible offline/online decomposition: First, a typically 
expensive offline-phase involving the computation of the reduced spaces is per- 
formed. This phase only needs to be precomputed once. Then, the online phase 
allows an extremely fast computation of the RB solutions for many new param- 
eters as only low dimensional systems need to be solved. Recen tly, we adopted 
the R B methodology to constrained stationary elliptic problems dHaasdonk et al , 



2011b . which we extend here to the instationary case. We refer to the recent contri- 



bution (ICont et all 1201 II) for a tailored RB approach in option pricing. In contrast 
to our setting no inequality constraints are taken into account. The main challenge 
is the construction of a suitable low dimensional approximation of the dual cone 
required for the approximation of the constraints. In this contribution, we introduce 
a new greedy strategy based on an angle criteria and show numerical results. 



2 American Option Model 

An American option is a contract which permits its owner to receive a certain 
payoff y(S, t) > at any time T between and T > 0. The variable T indi- 
cates the_jmturity^Inttodi^ng^ time variable t := T — X, we can use, 



e.g., dAchdou and Pironneau , 2005 ) the following non linear model 



d,P - ^a 2 s 2 d 2 s P - (r- q)sd s P + rP>0 7 P - y > 0, 
d t P - ^a 2 s 2 d; s P - (r- q)sd s P + rP^j ■ (P- yr) = 0, 

where P = P(s,t) is the price of an American put, with s G M + the asset's value, a is 
the volatility, r is the interest rate, q is the dividend payment and y = y/(s,t) is the 
payoff function. The boundary and initial conditions are as follows: P(s,0) = y(s), 
P(0,t) = K, lim. s ^+ooP(s,f) = 0, where K > is a fixed strike price that satisfies 
K = \if(Q,Q). In what follows, we focus on the case y/(s,t) = (K — s)+ with (•)+ = 
max(0, •), but our method applies as well to other types of payoff functions. For 
the implementation, we restrict the values of s to a bounded interval Q := (0,s/), 
where s/ is large enough to make the assumption P(sf,t) = realistic. Let us also set 
P = P Pq, with initial data Po(s,t) = K(l — s/sA, so that P satisfies homogeneous 
Dirichlet conditions. Our aim is now to reformulate the last system in a weak form, 
where our reduced basis method applies. In this view, we introduce the following 
functional spaces: 

V := {v e L 2 {Q)\sd s v £ L 2 {Q),v {dn = 0}, W := V' . 
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The scalar product (v)y associated with V is defined by (u,v)v'= ( s d s u,sd s v) L 2^ + 

(u,v) L 2(q\, where (■,-)l 2 (q) ^ s tne usual scalar product on L 2 (Q). The operators are 
specified as follows: 

1 , , 

a(u,v,H) = -a {d s u,d s (s v)) L i {n) + {-(r-q)sd s u + ru,v) L 2 {£2) , 

/(v;m) = (F^L^n), «(i?;m) = (9,*l)w, 

with F := -(^Po-^V^Po-^-^K^o + ^o), i-e. F = K(±q-r) and 
)// := y/— Pq. For rj € W = V', we also define b(rj,v) = f](v). We can now recast our 
problem in the following weak form, parametrized by n = (K,r,q,a) £ ^ C R 4 . 
We now introduce m as a weak representant of the solution P, as this is the standard 
notation in reduced basis literature: 

(d t u,v) L 2 {n) +a(u,v;^)-b(X,v) =/(v;/x), veV (1) 
2>(t]-A,w) > g(7j-A;M). neifef, (2) 

where M C W is a closed convex cone. Various methods can be considered to solve 
numerically Equations In what follows, we use a 0-scheme for the time dis- 

cretization. Given fi € L £ N and At := T /L, this method corresponds to the 
following iteration. 

Given < n < L - 1 and u" e V, find u" +l e V and A" +1 £ M that satisfy Vv £ 
V,V77 gM, 

/ U »+ -U n \ +fl(0M »+l + (1 _ 0)M » )V . M) _ fo(A n+l )V)=/(v;M)) (3) 

KT?-A n+1 , M " +1 ) >g(Tj-A" +1 ;M)-(4) 

This recursive definition is initialized with u° := Note that in this scheme, the 
definition of A" is not recursive. 



3 Reduced Basis Method 



Standard finite element approaches do not exploit the structure of the solution and 
for a given parameter value, a high dimensional system has to be solved. In what fol- 
lows, we introduce a specific Galerkin approximation of the solution, based on the 
reduced basis method and present algorithms to compute the corresponding bases. 
The principle of the reduced basis method consists in computing parametric solu- 
tions in low dimensional subspaces of V and W that are generated with particular so- 
lutions of our problem. Let us explain in more detail the corresponding formulation. 
For N £N, consider a finite subset := {jUi , . . . , Hn} C & with ^ fij, Vi ^ j. 
The reduced spaces Vn and Wn are defined by Vn := span { y\ , . . . , \j/^ v } and 
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Wn := span {^i , . . . , ^ w } where y/, and are defined from the large set of snap- 
shot solutions u n (ni) and A"(ji;), i = l,...,N, n = 0, . . . ,L. Here k"(jUi) and A"(jUi) 
denote the solution of Equations ([3}j4]i at the time t„ := nAt for the parameter value 
jj. = muj. The functions y/j and %j are suitably selected elements spanning Vjy and 
Wat with Ny,Nw <N(L+l) preferably small. Both families fV = (Wj) j=i,...,n v and 
SiV = (E,j)j=i,...,N w are supposed to be composed of linearly independent functions, 
hence are so called reduced bases. Numerical algorithms to build these two sets will 
be presented in Section|4] We define the reduced cone Mn C M as 

r % 

In this setting, the reduced problem reads: 

Given jtx G < « < L - 1, i$ G V N , find m^ +1 G Va< and A^ +1 £ M w that satisfy 
Vv N eV N ,Vri N eM N , 

\ At +a ^ +l + {l - 0)u n N ,v N ;n)-b(^ + \v N )=f(v N ;n), (5) 

\ / L 2 (£2) 

fo(^-A^ +1 ,4 +1 )>g(Tj w -^ +1 ;M),(6) 
where the initial value ajy is chosen as the orthogonal projection of uq on V^v, i.e. 

(u N - u , vtf) v = 0, Wiv G V N . 



4 Reduced Basis Construction 



In this section, we present two methods to extract a basis fV C V and En CM from 
the snapshots. Both are greedy procedures based on a finite training set £? t rain C 
small enough such that it can be scanned quickly. Given an arbitrary integer Nw, the 

dual reduced basis En = (^j)j=i n w is built iteratively according to the following 

algorithm. The goal of the approach is to obtain a reduced cone Mjy C M capturing 
as much "volume" as possible. 

Algorithm 1 (Angle- greedy algorithm) Given Nw, Strain C choose arbitrarily 
< n\ < L and Hi G Strain an d do 

1. set E N = { |ff ( [fij k }, W N : = span(S^), 

2. for k = 1 , . . . ,N\y — 1, do 

a. find (n k+1 ,flk+i) := argmax, 1=0 L ^,9> lmm U (^"0),W#)) , 

b seth^ — Vk+1 ^ 
D. setq k+1 . {{l n k+1(jlk+l)llw , 

c. define E k N +1 = S* U{^ +1 }, W N +l := span(4 +1 ), 
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3. define En '■= E^ w , Wn '■= span(SAr). 

Here we have used the notation /(v,5) to denote the angle between a vector v and a 
linear space S CW, which is simply obtained via the orthogonal projection lis from 
W on S by 



Z(v,5) = arccos ■ 



\ v \\w 



vew. 



We apply the POD-greedy algorithm dHaasdonk and Ohlbergerl 120081) to design 
the primal reduced basis fV This procedure is standard in RB-methods for evolu- 
tion problems. In RB-methods, frequently weak greedy procedures are used, which 
make benefi cial use of r a pidly computable error estimators and allow to handle large 
sets Strain ( IBuffa et all 1201 II) . However, as our analysis does not yet provide a- 



posteriori error estimators, we use the true projection errors as error indicators. This 
corre sponds to the so called strong greedy procedure ( IBuffa et all 1201 ll lHaasdonk , 
201 lb . 



Algorithm 2 (POD-greedy algorithm) Given Ny > 0, 2? train C choose arbi- 
trarily Hi e Strain, 

2. fork=l,...,N v -l, do 



a. define fl k+1 := argmax^,^ (£ B=0 \\u"(fl) - U^(u n (n))\ 

b. define ijf k+1 := POD \{u"{pL k+l )- JI^ (u n {pL k+ 1 ) ) 

c. define^ 1 :=^U{^ +1 }, 
3. define T N := , V N := span'fV. 

Here, we have denoted by Tlyk the orthogonal projection on with respect to 

N 

(■,-)v, an d by POD\ the routine that extracts from a family of vectors the first 
Proper Orthogonal Decomposition (POD) mode that can be obtained via the best 
approximation property 



POD i (v")„ = 



(i /. • 



argmin £ ||v* - (v",z) v z| 
M\v=i n=0 



2 

V 



In this definition V is spanned by v", n = 0,...,L. A con vergence analysis of the 



POD-greedy procedure is provided in dHaasdonkl 1201 ll) . Note that Algorithm [2] 
always returns an orthonormal basis. This is even the case if a parameter value 
H <G Strain is selected more than once. We point out that our System (IS}© has a 
saddle point structure. Thus taking spanfV as reduced basis for the primal vari- 
able might result in an ill posed problem. T o guarantee th e inf-sup stability of our 
approac h, we follow an idea in troduced in (IRozzal 120051) for the Stokes problem, 
see also ( lHaasdonk et a 1|) for variational inequalities. It consists in the enrich- 
ment W N := U (B^,)i=i n W ' where B & is the solution of b{^,v) = (B^, v)y, 
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for v G V. We conclude with the final reduced space Vjy := span^N of dimension 
Ny :=dim Vpf. By construction we have Ny < Ny < Ny +N\y- 



5 Numerical Results 

In this section, we present some numerical results obtained on the American Option 
model. We start with a description of the numerical values and methods we use. In 
order to compute snapshots, we use a standard finite element method for the space 
discretization and the 0-scheme presented in Section[2]for the time-discretization. 
The time domain [0, T) = [0,1] is discretized with a uniform mesh of step size At := 
T/L, L = 20. The 0-scheme is used with = 1/2, i.e. we apply a Crank-Nicolson 
method. The space domain Q = (0,Sf) = (0,300) is discretized with a uniform 
mesh of step size As := Sf/S, S = 101. For the function space, we use standard 
conforming nodal first order finite elements. For the sake of simplicity, we keep 
the notation V for the discrete high dimensional space and define it by V := {v G 
ffi(i2)|v| [sm >Sm+1 j G Pi , m = 0, . . . , S - 1 } of dimension Hy = H := S -2 = 99 with 
s,„ : = mAs. We associate the basis function 0,- G V with its Lagrange node s, G £2, 
i.e., 0j(s/) = 8jj,i,j = l,...,H. The discretization of the Lagrange multipliers is 
performed using a dual finite element basis %j of W := V' having the same support 
as (j>j, so that b((j)j,Xj) = ?Hj, i,j = L • • • >Hw = H. The cone M is defined by: M = 
{l^i ^iiXii Hi > o| . To build the basis, we consider a subset Strain of 3 s that is 
composed of N = 16 values chosen randomly in the set 

&=[(1- §)*>, (1 + f )Zb] x [(1 - f )ro, (1 + |)r ] 



§)9o,(l + f)9o]x[(l-f)ob,(l + §) 



with the numerical values e = 0.1, #b = 100, r = 0.05, q = 0.0015, CT = 0.5. 
To define the basis ¥V and the convex set En, we use Algorithm [2] combined with 
the enlargement by the supremizers and AlgorithmQ] The eight first vectors of fV, 
En and the supremizers are represented in Figure Q] We simulate two trajectories 
corresponding to the values (Ny,Nw) = (8,8) and (Ny,Nw) — (16, 16) respectively. 
The corresponding bases fV are of size Ny = 16 and Ny = 32 respectively. We 
chose randomly a parameter vector ji corresponding to the values K = 106.882366, 
r = 0.048470, d = 0.007679, a = 0.418561 in Some steps of the simulation 
are represented in Figure [2] the top and lower row refer to the smaller and larger 
reduced spaces, respectively. We clearly see the improvement in the approximation 
by increasing the reduced dimensions. In order to evaluate the efficiency of the 
greedy algorithms proposed in Section |4] we plot the evolution of the quantities 



L 



e N := max 



£||««(M)-n v ,(M»(M))||^4:= max U(X n (ji),w£ 
MG^™» V n^0 N n = 0,...,L, V V 

fl £ £P train 
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s i s 

Fig. 1 Eight first vectors of the reduced basis fV, Sjv and the corresponding supremizers. 



Time step =1 Time step =10 Time step =20 




s s s 



Fig. 2 Finite element approximation (solid red line) and Reduced basis approximation (blue +) at 
time steps t /At = 1, t /At = 10 and t / At = T / At = 20. The payoff function y is represented with 
the black dashed line. The reduced bases that are used have been generated by (Ny,Nw) = (8,8) 
(plots on the top) or (Ny ,Nw) = (16, 16) (plots on the bottom). 



during their iterations. The results are plotted in the first two diagrams in Figure [3] 
We observe an excellent exponential convergence of the approximation measures. As 
final experiment, we address the generalization ability of the RB-model to param- 
eters outside the training set. We consider ^test C s , a random set of N te st = 10 
parameter vectors and estimate, for a given jj. 6 3? ', the efficiency of our method 
through these quantities: 



err N (n) = * AtY\\u n (n)-u%(n)\\v, Errjf = max {err N (n)). 
V „=o 

Note that errN(jj.) actually depends on fV ; for the sake of simplicity, we have 
omitted this reliance in the notation. As a test, we evaluate the influence of the pa- 
rameters Ny,Nw determining the sizes of the bases fV and Sjv on Ertjj . The results 
are plotted in the right diagram of Figure [3] In our example we numerically obtain 
Ny = Ny + Nw in all cases, indicating, that the primal snapshots and supremizers 
are linearly independent. We observe a reasonable good error decay when simulta- 
neously increasing AV and Nw, indicating that the reduced method is working well. 
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We also note that in our case, the size of the dual basis has a limited impact on the 
results. 



Angle-Greedy POD-Greedy 

1 1 j 10 j 1 1 1 1 1 1 1 1 1 ] Max error 

io2 




Fig, 3 Values of e^J and during the iterations of the greedy Algorithms Q] (left) and[2](middle). 
Right: Values of Enj^ with respect to Ny and Nyr. 
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